%% Loop over different values of rho
clear 
clc


% usual setup

% personid quarter drugid chosen cost_sharing incumbent control_extra treat_extra
% starts in June 1996 (1)
% can go to May 2006 (40)
A = importdata('general_entry_choice_data.txt');

data = A(:,1:6);
id = A(:,1);
quarter = A(:,2);
drugid = A(:,3);
chosen = A(:,4);
cost_sharing = A(:,5);
incumbent = A(:,6);

% % filtering
filter = (quarter <= 40) .* (drugid < 4) .* (1-isnan(cost_sharing)); % for now, filter based on brand, pre-2006
id(filter == 0, :) = [];
quarter(filter == 0, :) = [];
drugid(filter == 0, :) = [];
chosen(filter == 0, :) = [];
cost_sharing(filter == 0, :) = [];
incumbent(filter == 0, :) = [];
data(filter == 0, :) = [];


% extra parameters: number of periods, etc.
num_drugs = max(drugid);
num_periods = max(quarter);
num_people = max(id);

% setup work common to each loop
[~,~,id] = unique([id quarter], 'rows'); % unique person/quarter combos => creates choice id
outside = 1 - accumarray(id, chosen); % chose the outside option


rho_grid = -0.5:0.25:0.5;

for k=1:length(rho_grid)

    % Test some code in terms of simulated MLE
    % person, period, drugid, choice, copay, prev

    rho = rho_grid(k);
    rho


    % draw some normal random numbers
    rng(17)
    
    % NEW: multivariate normal
    S = 100;
    %random_draws = normrnd(0,1,num_people*num_drugs, S);
    random_draws = zeros(num_people*num_drugs, S);
    
    mu = zeros(3,1);
    Sigma = [1 rho rho; rho 1 rho; rho rho 1];
        
    for jj=1:S

        R = mvnrnd(mu,Sigma,num_people);
        R_reshaped = reshape(R',[num_people*num_drugs 1]); % reshape reads by column, so transpose it first
        random_draws(:,S) = R_reshaped;
    end

    
    delta_init = [-1.0783;-0.3111;-1.6126];
    %params_init = [delta_init; 0.04; 5; 1];
    params_init = [delta_init; 0.026;5.239;-0.910];
    params = params_init;

    myopts = optimset('TolFun', 1e-10, 'MaxFunEvals', 12000, 'MaxIter',10000, 'display','iter'); % off, iter, final, notify
    % previously 10^-8

    % test
    SMLE_v2(params_init,data,random_draws,num_drugs,num_periods,id,outside,S);

    [params_opt, mle] = fminsearch (@(t)...
        -SMLE_v2(t,data,random_draws,num_drugs,num_periods,id,outside,S), params_init, myopts);


    delta = params_opt(1:3);

    alpha = params_opt(4);
    gamma = params_opt(5);
    sigma = params_opt(6);

    delta
    [alpha, gamma, sigma]



    %% Compute standard errors
    fun = @(t) -SMLE_v2(t,data,random_draws, num_drugs, num_periods, id, outside, S);


    % hessian of a function at the final point
    %[H_small, H, H_big, G] = own_hessian_v2(fun, params_opt);
    [H_small, H, H_big, G] = own_hessian_v3(fun, params_opt);

    n=length(params_opt);
    for i=1:n
        for j=i+1:n
            H(i,j) = H(j,i);
            H_small(i,j) = H_small(j,i);
            H_big(i,j) = H_big(j,i);
        end
    end

    %variance_mat = inv(H_big);
    variance_mat = inv(H_small);
    standard_errors = diag(variance_mat.^(.5));

    alpha_err = standard_errors(4);
    gamma_err = standard_errors(5);
    sigma_err = standard_errors(6);

    [alpha gamma sigma]
    [alpha_err gamma_err sigma_err]


    save(['smle_correlated_' num2str(k) '.mat']);
end
